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ABSTRACT 

We report on a numerical evaluation of the statistical bootstrap as a tech- 
nique for radio-interferometric imaging fidelity assessment. The development of a 
fidelity assessment technique is an important scientific prerequisite for automated 
pipeline reduction of data from modern radio interferometers. We evaluate the 
statistical performance of two bootstrap methods, the model-based bootstrap 
and subsample bootstrap, against a Monte Carlo parametric simulation, using 
interferometric polarization calibration and imaging as the representative prob- 
lem under study. We find both statistical resampling techniques to be viable 
approaches to radio-interferometric imaging fidelity assessment which merit fur- 
ther investigation. We also report on the development and implementation of a 
new self-calibration algorithm for radio-interferometric polarimetry which makes 
no approximations for the polarization source model. 

Subject headings: techniques: image processing — methods: statistical — tech- 
niques: interferometric — techniques: polarimetric 



1. Introduction 



Radio interferometry produces astronomical images of unmatched spatial resolution and 
of unique scientific value. The images are formed by solving an inverse imaging problem 
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connecting the spatial coherence of the incident radiation, gridded into the visibihty plane 
orthogonal to the source direction, and the radio brightness distribution projected on the 
plane of the sky. An extensive review of this observing method is provided by Thompson, 
Moran, & Swenson (2001). The techniques used for indirect imaging in this discipline are 
well-established and robust, when used within the realm in which they are mathematically 
applicable. Best practices have been established in the community for their practical appli- 
cation based on algorithm research and evaluation, and common data reduction experience 
(Perley, Schwab, & Bridle 1989; Taylor, CariUi, & Perley 1999). 

A central unsolved question, however, concerns the development of a quantitative method 
to assess the fidelity and precision of the reconstructed images at each pixel on the plane of the 
sky. Radio synthesis image formation requires a simultaneous solution for the source bright- 
ness distribution and instrumental calibration effects, including the effects of atmospheric 
propagation and the instrumental response of the antennas and receiving electronics of the 
system (Rogers ct al. 1974). This joint solution is generally solved iteratively (Readhead 
& Wilkinson 1978), with an image deconvolution step per iteration to remove the effects 
of the sparse sampling of the measured correlation data in the visibility plane (Hogbom 
1974). Errors in calibration and imaging both introduce errors into the final reconstructed 
images. Formally, the inverse imaging problem admits an infinite family of solutions for the 
source brightness distribution due to the sparse samphng of the visibihty data. However, the 
problem can be regularized by applying physically reasonable constraints, such as positivity 
and compact support, to provide a robust calibration and imaging method which is strongly 
convergent in practice (Thompson, Moran, & Swenson 2001). 

An analytic solution for the fidelity of the reconstructed image is not tractable given the 
complex instrumental model and the common use of non-linear deconvolution methods. A 
method for quantitative fidelity assessment is, however, of vital scientific importance in ra- 
dio interferometry. The increasing instrumental complexity and output data rates of current 
and future radio interferometers, such as the Atacama Large Millimeter Array (ALMA^), the 
pre-cursor Combined Array for Research in Millimeter-wave Astronomy (CARMA^), the Ex- 
panded VLA Project (EVLA^), and the Square Kilometer Array (SKA^), require automated 
pipelined data reduction, as opposed to custom interactive reduction, if these telescopes are 
to achieve their full scientific potential and be accessible to the broadest astronomical user 



^http://www.alma. nrao.edu 

^http: / / www.minarray.org 
■^http: / /www.aoc.nrao.cdu/cvla 
http://www.skatelescope.org 
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community. This is particularly important for potential observers from other wavebands, or 
those who do not have a significant investment in specialized radio interferometric expertise 
at their local institutions. Pipelined data reduction can remove these barriers to entry, but 
requires a method for quantitative fidelity assessment to allow the statistical significance of 
features in the final science images to be clearly delineated. This is also true for general sci- 
entific analysis and interpretation of archived radio data, in repositories such as the National 
Virtual Observatory (NVO^). Apart from their importance in pipeline and archive analysis, 
we also note that fidelity assessment methods are useful in standard interactive data reduc- 
tion, particularly for long-duration targeted imaging observations where high-fidelity is a 
specific goal. 

Recent advances in techniques in computational statistics and the availability of leading- 
edge community resources for high-performance computing (HPC), allow fundamentally new 
approaches to be taken to the problem of radio-interferometric imaging fidelity assessment. 
We report here on the investigation of one such modern statistical technique, namely boot- 
strap resampling of dependent data (Lahiri 2003), applied to the test case of polarimctric 
calibration and imaging of Very Long Baseline Interferometry (VLBI) data. This test prob- 
lem has been chosen for study as it is broadly representative of the calibration and imaging 
problems in radio interferometry in general, and is also of intrinsic scientific and technical 
interest. It is not, however, the only imaging application in radio interferometry for which 
fidelity assessment is scientifically important. Other such problems include, for example, 
high-fidelity imaging of strong, extended sources in total intensity, and we reserve their 
study for future work. We note the use of the bootstrap to assess the functional precision of 
medical tomography images in a related inverse imaging problem domain (Maitra 1997). 

The advent of the Very Long Basehne Array (VLBA^), which was engineered with a 
uniform instrumental polarization response and a higher overall sensitivity than previous 
US VLBI arrays (Kellermann Sz Thompson 1985), has significantly broadened the scope and 
range of science possible using VLBI polarimetry(Wardle & Roberts 1994). Polarization 
VLBI studies are frequently concerned with imaging weak polarized emission, as common 
target sources such as continuum extra-galactic radio source cores typically have a degree of 
integrated linear polarization not exceeding several percent of the total Stokes / brightness 
(Cawthorne et al. 1993). Image fidehty assessment methods are therefore particularly rele- 
vant to VLBI polarimetry (Roberts, Wardle, & Brown 1994). We note that although other 



^http://www.us- vo.org 

^The National Radio Astronomy Observatory is operated by Associated Universities, Inc., under cooper- 
ative agreement with the National Science Foundation 
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emission processes studied with VLBI, such as maser emission, may have a higher degree of 
either hnear or circular polarization than continuum synchrotron sources (Kemball 2002), it 
remains equally important to be able to assess the significance of weak polarized emission in 
imaging studies of these objects as well. 

In this paper we explore the use of the statistical technique of bootstrap resampling 
to estimate the pixel-based variance in the output images produced by the calibration and 
imaging of simulated polarization data from the VLBA. These variance images are inter- 
comparcd with results obtained by direct Monte Carlo simulation, and confirm the validity 
of the bootstrap approach in this domain. 

The paper is structured as follows: Section 2 describes the theory of polarization calibra- 
tion and bootstrap resampling. The details of the methods used in this numerical study are 
presented in Section 3 and the numerical results in Section 4. The discussion and conclusions 
are presented in Section 5 and Section 6 respectively. 

2. Theory 

2.1. The imaging equation for VLBI polarimetry 

Radio-interferometric polarimetry, which allows imaging of the radio brightness in all 
four Stokes parameters {I,Q,U,V}, is possible if a sufficient subset of antennas in the ar- 
ray is equipped with orthogonal polarization receptors (usually crossed-linear or oppositely- 
polarized circular) and the correlator forms all available polarization cross-products between 
the independent polarization channels at each antenna on each baseline. The calibration al- 
gebra for interferometric polarimetry was first derived by Morris, Radhakrishnan, & Seielstad 
(1964), with further work provided by Conway & Kronberg (1969) and Weiler (1973). The 
analysis for the specialization of continuum VLBI polarimetry has been developed systemat- 
ically over the past twenty years (Cotton et al. 1984; Roberts ct al. 1984; Roberts, Brown, & 
Wardle 1991; Cotton 1993; Roberts, Wardle, & Brown 1994; Leppanen, Zensus, & Diamond 
1995; Kemball, Diamond, & Pauhny-Toth 1996). A corresponding signal path analysis for 
spectral line polarization VLBI, including a formulation of the calibration and propagation 
effects along the signal path from the top of the atmosphere through the antennas, feeds, 
receiving electronics and correlator, was provided by Kemball (1993) and Kemball, Diamond, 
& Cotton (1995). 

A general mathematical framework for radio-interferometric polarimetry is described 
in a series of papers by Hamaker, Bregman, & Sault (1996), Sault, Hamaker, & Bregman 
(1996), Hamaker & Bregman (1996) and Hamaker (2000). This series of papers (hereinafter 
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HBS) parametrizes the signal path analysis at each antenna as the product of a series of 
(2 X 2) Jones matrices, by analogy with their use in transmission optics. This mathemat- 
ical formulation represents no new physics over earlier analyzes of the signal path and the 
resultant polarimetric interferometer response. However the use of Jones matrices and the 
outer matrix product allows a concise formulation in which the calibration algebra can be 
expressed independent of the choice of a specific polarization receptor basis, be that linear 
or circular. 

The outer matrix product for two Jones matrices A and B, denoted B, is cited by 
HBS as: 
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Each quadrant of A (E) B is formed as the product of the corresponding element in A 
times the matrix B. The outer product of two n x n matrices has dimension x n^. 

(Cornwell 1995a) generalized the HBS framework to include image-plane calibration 
effects, leading to an imaging equation of the general form: 



n ® ^n] I n ^T^^P) ® '^n{P)\ e-'^^'--^^-^-^^ K S{p)dn (2) 



where denotes the outer matrix product, j = y— T; S{p) is the radio brightness 
expressed as a Stokes four- vector toward the unit direction p & fl, and is the region of 
astronomical interest on the sky centered on pi, 
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Vmn is the measured visibility correlation between antennas in and n in a polarization 
receptor basis {p,q) G {{X,Y), {R, L)}, for linearly- and circularly-polarized feeds respec- 
tively. 
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is an element in the sub-series product of Jones matrices which have no direction 
dependence, ^ f{p)j thus describes visibihty-plane signal path corrections at an- 
tenna m of type k, where k enumerates the sub-types of the calibration corrections in the 
signal path. Examples of visibility-plane calibration correction sub-types enumerated by n 
include those arising from bandpass and electronic gain corrections. In the image plane, 

is an element in the product of Jones matrices with direction dependence, = f{p), 
thus describing image-plane corrections at antenna m of type k, and K is a fixed (4 x 4) 
conversion matrix which maps the Stokes four-vector S{p) into polarization correlation pairs 

in the receptor polarization basis. 

We omit decorrelation losses (HBS) and baseline-dependent digital signal processing 
corrections (Cornwell 1995a) in this expression, as they are not relevant to the work presented 
in this paper. The individual Jones matrix elements for each correction matrix are expressed 
in the adopted polarization receptor basis (e^, e^) and can also be parametrized arbitrarily 
on parameters defining an individual instrumental model: 



The sets of Jones matrices in the visibility-plane, {Gm}, and the image-plane, {Tm}, 
model the calibration corrections for the full signal path at each antenna. These individual 
calibration components are enumerated in further detail by Cornwell (1995a), Noordam 
(1995), and Noordam (1996). For the chosen case of VLBI polarimetry used here to evaluate 
the bootstrap resampling for radio-interferometric imaging fidelity assessment, we restrict the 
instrumental corrections to the Jones matrices for the parallactic angle correction, Pm = G^, 
and the instrumental polarization response. Dm = G^. These formally both depend on 
direction, p, but for the narrow field of view typical in VLBI we consider these as direction- 
independent corrections at the center of the field, with Pm,Dm G {Gm}- These are the 
on-axis values accordingly. For the VLBA antennas, which have altitude- azimuth (alt-az) 
mounts and circularly polarized receptors. 
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where am{t) is the time- variable parallactic angle at antenna m, which is known an- 
alytically (TMS). The instrumental polarization response, also known as the D-terms or 
polarization leakage (Morris, Radhakrishnan, & Seielstad 1964; Conway & Kronberg 1969), 
take the Jones matrix form: 



where dm are the traditional complex components of the D-term at antenna m for 
polarization {p,q) (Noordam 1995). 



By design of the presented work, the VLBl polarimetry data considered here require no 
calibration for electronic gain amplitude or phase. Equivalently stated, there is no diagonal 
Jones matrix in {Gm}, outside of or D^, which needs to be applied to correct the data. 
This is the usual condition for polarization VLBI data in current data reduction practice be- 
fore the final incremental solution for the polarization leakage terms. The residual problem of 
solving for Dm requires knowledge of S{p) or a joint solution for S{p) and Dm simultaneously, 
as is evident from equation 2. In general, this problem is more difficult for polarization VLBI 
than connected-element interferometry, due to the higher spatial resolution and the relative 
absence of calibrator sources with a known polarization structure at milliarcsccond (mas) 
resolution. A number of approximations for the source model polarization structure have 
been developed in the past to sharply reduce the number of free parameters x/ describing the 
source polarization model, and thus allow a single fit for the source model unknowns xi and 
the Dm simultaneously. For an alt-az array such as the VLBA, each antenna has a sufficiently 
different variation of parallactic angle with time to allow these parameters to be separated. 
Expressed equivalently, the basis functions e-'"'"^*^ are non-degenerate over a sufficient range 
in parallactic angle variation. The source approximation methods in common VLBI use 
are summarized by Kemball (1999). These approximations include assumptions of: i) a 
linearly-unpolarized calibrator, {Q,U) = (Roberts et al. 1984); ii) an unresolved calibra- 
tor, (Q, U) = >c (Cotton et al. 1984); iii) a similarity approximation, Q + jU = xJ (Cotton 
1993); iv) a multi-component similarity approximation, ^i{Qi + jUi) = xili (Leppanen, 
Zensus, & Diamond 1995), and, v) a spectral- line approximation, Q{i') -\- jU{v) — x^/(i/) 
(Kemball and Diamond 1997). This family of solution methods will be referred to as source 
approximation methods in what follows. These approximations reduce the source polariza- 
tion model parameters to fewer than required for a fully general polarization model by fixing 




(7) 



2.2. Strategies for VLBI instrumental polarization calibration 
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the scaling factor x per image (iii), per Gaussian image component (iv), or per spectral 
channel (v) respectively. 

It has also been common practice in polarization VLBI to linearize the feed leakage 
cahbration equations by ignoring terms of order O(D^) or 0{D • {Q + jU)) (Roberts et 
al. 1984) and to use only the resulting hnear equations for the cross-polarized visibilities 

(ypq^ yqp^ in a linear fit for and D,n (Cotton 1993; Kcmball, Diamond, & Cotton 1995; 
Leppanen, Zensus, & Diamond 1995; Kemball and Diamond 1997). This linearization has 
also typically been applied to the equation inverse when correcting the visibility data for 
the polarization leakage D— terms. The use of the matrix-based HBS framework implicitly 
includes all terms without linearization. 



2.3. Polarization self-calibration 

Polarization self-calibration has been proposed frequently in the literature as a likely 
optimal technique for solving simultaneously for the source model and the instrumental po- 
larization response (Roberts, Wardle, & Brown 1994; Sault, Hamaker, & Bregman 1996). 
This is by direct analogy with the cornerstone role self-calibration plays in radio interfer- 
ometry in allowing a simultaneous solution for the total intensity distribution and complex 
antenna-based electronic gains. A recent review of electronic gain self-calibration is provided 
by Cornwell & Fomalont (1999). Several of the source approximation methods described in 
the previous section have already been applied iteratively, as a variant of polarization self- 
calibration (Kemball, Diamond, & Pauliny-Toth 1996; Kemball and Diamond 1997). 

A simple definition of polarization self-calibration is that of an joint (often iterative) 
solution for S{p) and within a calibration framework for radio-interferometric polarime- 
try, with optional constraints applied in the visibility and image planes to regularize the 
solution as required. As such, it is the direct analog of standard total intensity phase and 
amplitude self-calibration. In the HBS framework, general self-calibration is simply any joint 
solution for a set of unsolved Jones calibration matrices in {Gm} and {T^} simultaneously 
with the source brightness distribution S{p)] all are free parameters in the imaging equation 
on a equal footing. In this nomenclature, polarization self-calibration, as defined above, is a 
simple subset of general self-calibration in the imaging equation formalism. 

In the imaging equation framework [2] , given an estimate for and the known analytic 
Pmi the source brightness can be determined from and the use of standard deconvolu- 
tion techniques, as described by (Cornwell 1995b), where is computed as the difference 
between the observed and predicted visibility data. Conversely, given an estimate of S{p), 
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and the known analytic Pm, can be formed at the position of the unknown Dm in the 
imaging equation by computing contributions from the right-hand and left-hand sides of [2] 
respectively through to the position of Dm in Hl^m ^ ^nil forming by inserting the 
current value of the D^ matrix at that position (Cornwell & Wieringa 1996). This also allows 
a direct computation of which can be used in a non-linear fit for the off-diagonal Jones 
matrix elements of D^. in each polarization solution interval. We re-express this mathemat- 
ically by defining two operators, TZ^ and which operate from the right- and left-hand 
sides of the imaging equation respectively, through to the position of a visibility calibration 
component being solved for. The operators yield 4-vectors containing the accumulated 
product of visibility cross-correlations at that index point in the imaging equation. For 
baseline mn, 
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The calibration solver minimizes. 



mn 

Polarization self-calibration can be achieved in this framework using an iterative solver 
which starts with an initial value, D^^, for each D Jones matrix, and iterates with optional 
visibility- or image-plane constraints to solve successively for S{p) and the set of Dm per 
antenna and polarization solution interval, within the imaging equation model. The data 
are pre-averaged in the fit over time intervals short with respect to the maxima | ^^^^^ | and 
l^^^^l where t is used here as a parameter along the uv-tracks in the visibility plane. We note 
that is formed here from the visibility 4-vectors, including cross-polarized correlations. 
The imaging deconvolution similarly uses the 4- vector of full polarization intensities. 



2.4. Imaging fidelity assessment 

The visibility data measured at loci in the uv-plane for each projected baseline mn 
constitute a time-series of complex observed data, denoted here as: {V^^{ukiVk,tk), k — 
1, Nmn, V (m, n)}. We adopt a simple additive noise model for the measured visibility data: 



(11) 



where V„in is given by equation [2], and A/" is a random phasor drawn from the complex 
normal probability density for CJ\f{0, of zero mean and variance cr^^. We adopt the 
statistical nomenclature here of Kay (1993) for complex random variables. The noise contri- 
bution is independent and identically distributed (HD). For a radio interferometer comprised 
of uniform array elements of system equivalent flux density SEFD (a fundamental parameter 
of antenna sensitivity), bandwidth Az^, and a sample integration time At, the thermal noise 
variance is given in the form (TMS): 



= 2A^ ^^^^ 

The complete set of measured visibihty data {V^^{uk, Vk, tk), k — 1, Nj^n-, V (m, n)} can 
be regarded as a single statistical realization of a sequence of random variables {Vi, V2, V^v} 
of sample size N from a joint multivariate probability density function (PDF), 



where the covariance matrix Cy = (^th^^ superscript ^ denotes conjugate transpose (Kay 
1993), U is an {N x TV) unit diagonal matrix, and the imaging equation [2] yields Vmn — 
f{S{p), C^, T^). Although this is a parametric distribution, here in terms of {S{p), G^, T^), 
these quantities are not known at the time of observation, and are determined by the calibra- 
tion and imaging process used in data reduction. We consider any solvers for the calibration 
matrices (G^jT^) or the source brightness distribution S{p) as statistical point estimators 
for these parameters and denote these estimators as S{p), G^, and T^. This perspective con- 
nects the observed visibility data {V^^{uk, Vk,tk), k = 1, N^nn, V (m, n)} and the unknown 
source brightness distribution and calibration parameters appearing in the PDF, in the form 
of a standard statistical inference problem. In this context, the problem of imaging fidelity 
assessment is equivalent to that of determining the sampling distribution properties of the 
estimators S{p),G'^, and T^. The imaging problem is assumed regularized, and therefore 
convergent and robust, in this analysis, in order to remove the poor conditioning implicit in 
the sparse sampling in the uv-plane. We omit the explicit choice of basis representation for 
S{p); this could be pixel-based or of functional form. 

We denote the sampling distribution of S{p) as ^s{p)- parent distribution 

p{V; SlpjjC^jT^) for the observed visibility data is unknown, the sampling distribution 
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J's{f) is also unknown a priori. We are physically and logistically constrained from measuring 
p(V; S{p), C^jT^) by drawing all potential realizations of the visibility time series from the 
parent population; each observing run is unique in terms of the instrumental and atmospheric 
conditions, quite apart from the logistical impracticality of this approach. In addition, if the 
calibration corrections are imperfectly modeled in any significant sense, = f{^i), the 
parent PDF will differ in mathematical form from the assumed model p{V; S{p),G'^,T^). 
Errors in the form of the parametric model compromise the statistical inference based on 
that model. 



2.5. Bootstrap resampling 

Bootstrap resampling techniques offer an alternative statistical inference method which 
is not as sensitive to errors in the assumed model and does not require that the problem 
be analytically tractable. This method constructs statistical resamples from the measured 
data realization which mirror the statistical properties of the unknown parent population. 
Developed by Efron (1979) for IID data, this method has proved to be a powerful technique 
in modern computational statistics, with increasingly broad applicability. It has undergone 
significant theoretical development since its inception and modern reviews are provided by 
Shao & Tu (1995), Davison & Hinklcy (1997), Chernick (1999), and Lahiri (2003). Bootstrap 
techniques have also been extended to the case of statistically dependent data which are not 
IID. This question is considered in detail in a monograph by Lahiri (2003), and is of particular 
importance for our study. The observed visibility data in radio interferometry are generated 
by a process with long-range statistical dependence, defined to first order by the Fourier 
transform of the source brightness distribution to the visibility ■uv- plane [2] . Other bootstrap 
resampling techniques, such as the jackknife and general subsampling methods (Shao & Tu 
1995; Politis, Romano, & Wolf 1999), are also relevant to resamphng of dependent data 
(Lahiri 2003). 

A single measured realization {V^^{uk,Vk,tk), k — 1,7V^„, V (m, n)} of the random 
process can be used to construct a distribution function ^(V; 5'(p), G^, T^) from which sta- 
tistical resamples V* can be drawn. The bootstrap principle requires that the statistical re- 
lationship between V* and p(V; S{p), C^, T^) mirror that between V and p(V; -S'(p), G^, T^) 
(Lahiri 2003). We consider the specific choice for ^(V; S{p), G^, T^) in further detail below. 

Wc denote resamples drawn from p{V; S'(p),G^,T^) conditional on the observed data 
{V:^':{uk,Vk,h), k = l,Nrnn, V (m,n)}^ as {V;^'^{uk,Vk,h), k = l,iV„„„ V (m,n)}t..^j 
when acted on by the imaging estimator S{p), they yield S^f. The bootstrap estimate of the 
sampling distribution, ^g/g) can then be obtained from the distribution of S^f. 
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The primary statistical requirement on successful use of the bootstrap method is that 
the relationship between the resamples and the PDF derived from the data realization accu- 
rately reflect the corresponding relationship between the parent population and the observed 
data. For IID data, the empirical probability distribution function is an optimal choice: 
p{») = X^fcLi II ("1^ ^ •)) where the indicator function \\{A) = 1 if A is true, else zero 
(Lahiri 2003). For IID data, resamples can be drawn from the empirical distribution func- 
tion directly, with replacement, to generate a valid bootstrap ensemble (Davison & Hinkley 
1997). 

The IID bootstrap is not statistically valid, however, for the case of dependent data 
Singh (1981), and modifications to address dependence are reviewed by (Lahiri 2003). We 
summarize them here in several broad categories: (i) block bootstraps; (ii) model-based 
bootstraps, and; (iii) subsampling bootstraps. Block bootstrap methods resample the data in 
blocks of sufficient length to capture the bulk of the auto-covariance dependence information 
in the data. Model-based bootstrap methods model and remove the statistical dependence 
from the data and apply the direct bootstrap to the residuals, assumed then to be IID 
(Lahiri 2003). In this bootstrap method, the model fit is tailored to the process generating 
the data. Subsampling draws unmodified fractional resamples from the data; as such it 
is a generalization of the delete-d jackknife (Politis, Romano, & Wolf 1999). Subsampling 
preserves directly the long-range statistical dependence in the resamples. 

Although the bootstrap is a computationally intensive technique, geometric advances 
in available HPC resources make this technique readily applicable as an approach to ra- 
dio interferometry imaging fidelity assessment. We discuss the practical evaluation of this 
technique with radio interferometry data in subsequent sections. 

3. Simulation methods 

We have conducted several numerical studies to inter-compare the results of bootstrap 
resampling with those obtained by direct Monte Carlo simulation, as a means to evaluate 
the applicability and usefulness of bootstrap resampling techniques in radio interferometric 
imaging fidelity assessment. As noted in the introduction, we have chosen VLBI polarization 
self-calibration and imaging as the test case for evaluation, both as it is broadly represen- 
tative of general calibration and imaging in radio interferometry, and also because of its 
intrinsic scientific and technical interest. In the language and formulation of Section 2.4, we 
seek to measiuc the sampling distribution properties of the imaging estimator using several 
candidate bootstrap resampling techniques, in order to assess their statistical performance 
in this problem domain. 
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3.1. 



Polarization self-calibration heuristic 



For the numerical studies conducted here, we have developed a polarization self-calibration 
algorithm based on the general principles outlined in Section 2.3. We do not consider the de- 
tailed scientific optimization of the polarization self-calibration algorithm here, as this aspect 
is not central to our study of the bootstrap method. We require only that the algorithm be 
broadly representative of typical calibration and imaging data reduction processes in radio 
interferometry. 

We adopt a simple polarization self-calibration heuristic accordingly, starting with an 
initial unit diagonal D Jones matrix at each antenna. 



and iterating without any visibility- or image-plane constraints to solve successively for 
S{p) and each per antenna, varying only the off-diagonal D Jones terms. The data 
were pre-averaged to 5 seconds by the solver, and the D Jones solutions were assumed to 
have no time-dependence. CLEAN deconvolution was used during imaging, with a stopping 
threshold equal to the expected thermal noise limit; this deconvolution threshold was held 
constant for all self-calibration cycles. In this implementation, ten iterations of calibration 
and imaging were performed using non-progressive self-calibration, in which the uncorrected 
observed data were used in each cycle, as opposed to the corrected data from the previous 
cycle. 

Iterative polarization self-calibration has not been used extensively in VLBI polarimetry 

partly due to the lack of software support in existing data reduction packages. The basic 
AIPS-I--I-'' package implements the imaging equation [2] and the HBS framework. For the 
applications developed here, we modified a reference version of the AIPS++ code base*^, held 
constant for reproducibility across the course of this work, and used the calibration solver 
and imaging implementation available in the package at that time. Additional utilities for 
bootstrap resamphng and analysis were also developed. These capabihtes are intended for 
future public release in a local package. 



''http://aips2.nrao.edu 
«vl.8 # 667 from 2003 
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3.2. Fidelity simulations 

3.2.1. Simulation configuration 

For these studies we chose to simulate single-channel data from the VLB A, at an observ- 
ing frequency of 43.0 GHz. Data for a source at (0:2000 — 12'*, (^2000 — +80^^^) were generated 
for all ten VLBA antennas in a set of 15 min scans starting at {tk — {k — 1)'* UT, k — 1, 20}, 
on an adopted date of 1 Feb 2003, using a correlator integration time At = 2 sec. The 
visibility data were generated in full polarization in a circular basis {RR, LR, RL, LL}. in- 
corporating arbitrarily chosen instrumental D-terms, which are tabulated for each antenna 
in table 1. The D-terms were chosen randomly but have a similar mean order of magnitude 
as those expected for the VLBA as a whole at this frequency. The source brightness distri- 
bution was chosen randomly as a pair of unit flux density elliptical Gaussian components 
matched to the resolution of the array, with properties summarized in table 2, and plotted in 
Stokes {/, Q, U, V} in figure 1. Their integrated polarization intensities in each polarization 
were also chosen at random, but arc of the same order as might be observed for 43 GHz SiO 
maser components using the VLBA. The parallactic angle of the Pie Town (FT) antenna at 
the center of the array varied from -21°, through -176° near transit, to -|-33° at the end of 
the simulation. This is plotted in figure 2. The array resolution for this synthesis was ~ 150 
Has in uniform weighting. Additive thermal noise was computed for each integration interval 
[11], using a nominal SEFD of 1436 Jy pubhshed for the VLBA 7mm band^. The data were 
generated at an image signal-to-noise ratio (SNR) of approximately 1300:1, as defined by 
the ratio of the peak source brightness to the off-source root-mean-square (rms) noise. This 
sensitivity is equivalent to that obtained by observing the unit fiux density source compo- 
nents in an 8 MHz bandwidth, or the same source components increased in flux density by 
a factor of 16, in a typical SiO channel bandwidth of ^ MHz. We adopted a normalized 
flux density scale for the source model in this work to allow different numerical results to be 
inter-compared readily, and a direct calculation of the dynamic range as the inverse rms. 



3.2.2. Monte Carlo simulation 

As the flrst step in the evaluation of bootstrap resampling, a reference statistical sam- 
ple containing A^^ = 256 realizations was generated by Monte Carlo simulation using the 
parametric model [13] and the array simulation configuration described above. The re- 
sulting Monte Carlo sample consisted of 256 simulated observed visibility time series data 



http: / /www. aoc.nrao.edu/vlba/ obstatus/obssum/obssum.html 
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{V^^{uk,Vk,tk), k — l^Nmn: V {m,n)}f^%^. The data were generated using a modified 
version of tlie task DTSIM, originally developed by the author in the AIPS^° package, so 
that the data would be generated and reduced in separate packages. Each realization was 
reduced using the polarization self-calibration and imaging heuristic described above, which 
is equivalent to the application of the imaging estimator S{p) defined in Section 2.4. We 
chose S{p) to represent the final restored image in a pixel basis Sxy, and computed estimator 
statistics for S{p) in this basis over the set of Ng restored images obtained from apphcation 
of the imaging estimator to the full Monte Carlo sample. A fixed circular restoring beam of 
full-width half-maximum (FWHM) equal to 156.007 fxas was used for all images, set to the 
geometric mean of the beam in uniform weighting for the ■uv-coverage of the simulated data, 
which is shown in figure 3. All images were formed of size (256 x 256) pixels with a pixel 
spacing of 30 ^as. 

The imaging estimator bias and mean-squared error (MSE) were computed per pixel 
with respect to the known source brightness model convolved with the fixed restoring beam, 
as shown in figure 1. The mean and variance of the sampling distribution ^s(p) 
estimator S{p), were similarly computed per pixel over the sample of A^^ restored images. 
The variance was computed using the approximation A^^^ '^{S'^y) ~ Sly, for computational 
efficiency. The MSE per pixel was computed as A^^""-*^ ^ {S^y — S^y) , where S^y is the known 
true source brightness model convolved with the fixed restoring beam. Each statistic was 
computed per Stokes polarization {/, Q, U, V} and per polarization self-calibration iteration 
number, across all Ng realizations. Similar estimator statistics were also accumulated per 
antenna and self-calibration iteration number for the polarization calibration estimators 
for each D Jones matrix determined by the polarization self-calibration and imaging solvers. 

The estimator sampling distribution properties estimated using Monte Carlo simulation 
were used as the reference in assessing the statistical performance of the bootstrap techniques. 
In practical radio interferometry, direct Monte Carlo simulation of this form is not practical, 
as there is no a priori knowledge of the exact source brightness distribution and calibration 
matrix values. 



3.2.3. Bootstrap simulations 



[MC 
1127 

from the Monte Carlo sample as input to the bootstrap resampling numerical studies. We 



We randomly chose the 127*'* realization {V^^{uk,Vk,tk), k = l,iV^„, V {Tn,n)}{2-j 



' http : / / www. aoc.nrao.edu /aips 
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resampled from this template realization using both model-based bootstrap methods and 

subsampling techniques, and generated in each case a sample of size A^^ = 256 resampled 
observed visibility time series {V^^{uk,Vk,tk), k = 1, A'"^^, V {'>TT-,n)}l The block boot- 
strap was not evaluated in this study due to the known long-range dependence in the data. 
The polarization self-calibration imaging estimator S{p) was applied to each realization and 
pixel-based estimator statistics computed across the S*y in the same manner as were com- 
puted for the Monte Carlo sample described above. These statistics estimate the bootstrap 
sampling distribution ^s(p) ^® inter-compared with those obtained from the Monte 

Carlo sample to evaluate the statistical performance of the bootstrap. 

For the model-based bootstrap, we adopted a piece-wise polynomial model for the real 
and imaginary parts of the observed visibility time series as a model of the long-term depen- 
dence in the data. A model of the form: 



\ 1=0 1=0 / 

was fit over successive bootstrap model intervals Atb on each baseline, using least- 
squares minimization and a conjugate- gradient solver as implemented in the OptSolve-|--|- 
package^^. The lowest polynomial degree was used in each solution interval which yielded 
convergence. The model residuals. 



. 1=0 1=0 



were re-centered to avoid introducing bias (Lahiri 2003; Davison & Hinkley 1997), 



and resampled, uniformly and randomly, with replacement, to yield V^^'^'^" {uk,V}.,tk)- 
The bootstrap resample for each /\ti, was then constituted as: 



V:^^{uk, Vk, tk)=[Y^ aiti + J Yl biA + VZ''''"iuk, Vk, tk) (18) 



, 1=0 1=0 



OptSolve-|— I- is distributed by Tech-X corporation (http://www.techxhome.com) 
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Four model-based bootstrap runs were performed, which are labeled Ml through M4. 
The bootstrap parameters used for each run are summarized in table 3. 

For the subsampling method, visibility points were deleted from the template realization 
randomly, avoiding repeated points, until a specified fraction, fg, of the total number of 
visibility points in each sample were removed for each subsample realization generated. Each 
subsample realization was processed using the imaging estimator S{p) as used previously, 

and the same estimator statistics accumulated over the sample of restored images as were 
computed for the model-based bootstrap and Monte Carlo runs. Four subsample bootstrap 
runs were performed, labeled SI through S4. The parameters used for each subsample run 
are summarized in table 4. 



3.2.4- HPC implementation 

All of the preceding Monte Carlo, bootstrap and subsample imaging estimator runs were 
computed on the public HPC Linux clusters deployed as a community resource by NCSA. 
Many of the bootstrap sample generation runs were also run in an HPC environment, but 
some were executed on a single workstation where possible. The parallel imaging estimator 
runs were mapped across the cluster nodes by assigning subsets of realizations to client nodes, 
and combining the partial statistical accumulations from each client node serially on a single 
node at the end of the run. A large degree of parallelism was achieved in these runs due to 
the limited communication needs between client nodes in this application. We typically ran 
over 32-64 cluster nodes, in individual runs of 4-8 hours duration. The duration depended 
on the underlying single-CPU performance of the specified cluster used. We expect these 
runs to scale efficiently to larger number of cluster nodes as they are well-matched to a 
loosely-coupled HPC architecture of this type. 

4. Simulation results 

4.1. Monte Carlo simulations 

For the Monte Carlo simulations discussed in the previous section, we plot the total 
MSE for the calibration estimator G^, summed over all antennas m, against polarization 
self-calibration iteration number in figure 4. 

The MSE, bias, and variance of the pixel-based imaging estimator S{p) are plotted for 
the restored Stokes {Q,U} images, at self-calibration iteration numbers 1 and 10, in figure 5 
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and figure 6 respectively. 



4.2. Bootstrap simulations 



To assess the statistical performance of the model-based and subsampling bootstrap 
methods, we have computed a goodness-of-fit statistic which compares the variance image, 
var^y, obtained by the imaging estimator for a particular bootstrap method against the 
corresponding reference variance image derived from Monte Carlo simulation, var^'-^ . For 
the subsample bootstrap, the estimator statistics need to be scaled by a factor expected 
to be proportional to ~ N{1 — fg)"- (Davison & Hinkley 1997), as derived for the delete-d 
jackknife, where the exponent, a, is dependent on the details of the individual model. We 
estimated the mean scaling factor, vj^ by summing over the inner quarter and over all Stokes 
polarizations {/, Q, V} of the ratio image , 



A goodness-of-fit statistic, vmse, was then computed over the same image region and 
polarization set, as: 



These values are summarized for the model-based and subsample bootstrap runs in ta- 
ble 5. For reference, the bootstrap parameters for the model-based and subsample bootstrap 
runs are summarized in table 3 and table 4 respectively. 

The variance images in Stokes {Q, U} for the final polarization self-cahbration iteration, 
obtained from the bootstrap runs, are plotted in figure 7, figure 8, figure 9, and figure 10. 

The variation of scaling factor, Vf, with subsample delete fraction, fg, is plotted in 
figure 11. 




(19) 




(20) 
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5. Discussion 

The polarization self-calibration algorithm developed as part of this study of the boot- 
strap method, has a convergence rate for these data, shown in figure 4, that is steeper than 
exponential as a function of iteration number. The scope of the current study does not allow 
us to extrapolate the quantitative statistical properties of this calibration estimator to all 
potential practical applications, but we believe this to be a viable, general technique for solv- 
ing for instrumental polarization in radio interferometry. Most importantly, it requires no 
polarization source model approximation, with the attendant systematic errors which may 
arise in this case. As is common for all polarization calibration algorithms in this class, we 
expect the solver to be sensitive to the degree of mutual degeneracy of the parallactic angle 
basis functions at each antenna (Conway & Kronberg 1969; Leppanen, Zensus, & Diamond 
1995; Sault, Hamaker, & Bregman 1996), which can invariably be equivalently quantified 
in terms of the maximum range in parallactic angle coverage for the array as a whole. As 
noted earlier, the simulated data used in this study were chosen to have good coverage in 
parallactic angle. 

Intuitively, we expect rapid convergence for polarization self-calibration algorithms from 
considering simple information arguments. For an interferometric observation of Nai self- 
calibration intervals, using an array of Nant elements, the number of self-calibration unknowns 
is N^tNant, and the number of knowns, defined by the data, is ■^«-^°"*(^^""*~^) ^ where Nt is 
the number of unique time-stamps in the observed visibility data. The ratio of knowns to 
unknowns is therefore r = ^*^2N^~^^ • -^^^ conventional amplitude and phase self-calibration, 
Nt ~ A^Ai, and r ~ ^^y^. For polarization self-calibration, the D-terms have slowly- varying 
time-dependence, A^a* ~ 1) and r ~ ^tNant ^ "Yhis ratio is significantly higher than for 
amplitude and phase self-calibration and we expect a steeper rate of convergence as a result. 
As a corollary, given that the number of free parameters being solved for in polarization 
self-cahbration is small relative to the number of visibility data points, it is more difficult 
to absorb coherent source brightness distribution errors in the calibration corrections than 
in conventional amplitude and phase self-calibration. This is dependent on the range of 
parallactic angle coverage however. 

This paper has considered the common case of constant, visibility-plane polarization 
calibration, but this polarization self-calibration algorithm can, in principle, be used to solve 
for time- and direction-dependent instrumental polarization corrections. This is required 
for high-fidelity, wide-field imaging studies but is not a routine mode of radio interferometer 
calibration at present. This is a subject for future study but we anticipate that the algorithm 
performance in this case will depend on the accuracy and number of free parameters in the 
parametrization adopted for the time- and direction-dependence of the D Jones polarization 



-20- 



corrections. 

Further information on the properties of the polarization self-cahbration algorithm can 
be obtained by examining the imaging estimator statistics for the Monte Carlo sample, 
plotted in figure 5 and figure 6. At the first iteration, D Jones matrix errors dominate, and 
the magnitude of the estimator bias is higher across the image as a result. At this iteration, 
the imaging estimator MSE is dominated by the estimator bias. At the final self-calibration 
iteration the D Jones matrix corrections are better constrained, and the bias is accordingly 
reduced. At this point the imaging estimator MSE is dominated by the estimator variance. 
Both calibration and deconvolution errors contribute to the measured imaging estimator 
statistics. To first order, errors in the estimated D Jones calibration corrections corrupt the 
cross-polarized visibility data with proportional fractions of the Stokes / visibility (Roberts, 
Wardle, & Brown 1994), so the variance in Stokes Q and U does not scale linearly with the 
underlying Stokes Q or U intensity at each pixel. 

The statistical performance of the model-based bootstrap, over runs Ml to M4, can be 
assessed by examining the goodness-of-fit statistic, vmse, tabulated in table 5, as well as the 
imaging estimator statistics plotted in figure 7 and figure 8. We expect optimal performance 
for the model-based bootstrap when the bootstrap model interval, Atb, is sufficiently long for 

adequate SNR in the fit to the bootstrap model parameters, but not so long as to invalidate 
the functional model for the statistical dependence in the data. The goodness-of-fit statistics 
listed in table 5, which capture the degree to which the bootstrap imaging variance matches 
that obtained from the reference Monte Carlo sample, suggest that run M3 reflects the 
optimal choice of model-based bootstrap run parameters for the data used in this study. 
Run M3 has the lowest vmse] therefore its variance estimates have the best overall fit to the 
image variance obtained by direct Monte Carlo simulation. 

It is clear, however, both from the tabulated goodness-of-fit statistic Vmse, and the 
imaging estimator statistics, that the model-based bootstrap performs well over a fairly broad 
range of Atb and Np] the values of these parameters used for the runs Ml to M4 are tabulated 
in table 3. We find this technique to be relatively robust with respect to the detailed 
bootstrap model assumptions, and therefore easy to use in practice. This robustness over 
a wide range of tuning parameters is a significant advantage of the model-based bootstrap. 
There are, however, differences between the model-based bootstrap resamples. Ml to M4. 
As this is a numerical simulation study, we know the parent population distribution [13], 
and can compare the empirical distribution of the resampled data against the known parent 
distribution. The closer the match the better the performance of the underlying bootstrap 
method. We plot this in the form of a cumulative distribution function (CDF) comparison 
in figure 12, for the real and imaginary parts of the data separately, for a randomly chosen 
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visibility row number 1000, in polarization cross-correlation LR. The empirical CDF for 
the reference Monte Carlo sample is plotted in figure 13. We denote the empirical CDF for 
the resampled data as c*^^'^^(V), and the true parent CDF as c^'^'^^(V), after removal of the 
known value of at this visibility point, and where, 

c^^''Hv(^'^)) = ^(l + erf(^— )) (21) 
^ V 2 ath 

The real and imaginary parts of the visibility data are denoted by superscript (i?, /) 
respectively. 

A goodness-of-fit statistic can be computed in the form cj^^^ = (c^^'^^ — c^^'^^y , 

and we tabulate tlic geometric mean of the value for the real and imaginary parts of the data 
cmse — s/cmse ^mse table 6. This analysis demonstrates that the bootstrap resamples 
generated for run M3 match most closely the parent distribution, consistent with our earlier 
result concerning the relative performance of the model-based bootstraps used in the current 
study. 

We note that the model-based bootstrap does require a detection of the bootstrap 
model in the basclinc-bascd bootstrap model intervals /\tb- This is broadly equivalent to the 
requirement that the data have sufficient SNR to allow conventional amplitude and phase 
self-calibration over a comparable interval. The bootstrap model is local to each interval At^ 
however, and is not equivalent to an overall model-fit to the source brightness distribution 
and cahbration corrections for the data as a whole. 

The statistical performance of the subsample bootstrap runs, SI to S4, can be de- 
termined by examination of the goodncss-of-fit statistic Vmse is table 5, and the imaging 
estimator statistics plotted in figure 9 and figure 10. From these results it is clear that the 
statistical performance of this bootstrap is more sensitive to the bootstrap run parameters, 
here defined by the subsample delete fraction fg, than the model-based bootstrap. Improved 
statistical performance is obtained for /s > | for the data used in this numerical study. This 
is consistent with the expected theoretical properties of the subsample bootstrap, which re- 
quires that — i> 1 as the sample size N ^ oo, if the bootstrap is to remain valid (Davison 
& Hinkley 1997). 

The variance scaling factor Vf is plotted against subsample delete fraction fs in figure 11. 
From theoretical considerations for the delete-d jacknife, we expect Vf oc (1 — /s)~" (Davison 
& Hinkley 1997). For the data used in this numerical study, a relation Vf oc (1 — fs)~^ agrees 
with the data to first-order, consistent with the expected theoretical result. 

Although the subsample bootstrap appears to under-perform the model-based boot- 
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strap, as assessed by the goodness-of-fit statistic vmse in table 5, this bootstrap has the 
advantage of requiring no fit to a local baseline-based bootstrap model and makes no as- 
sumption about the functional form of the long-range statistical dependence in the data. As 
such, it is a versatile technique which can find practical application in radio-interferometric 
imaging fidelity assessment. 

The computing costs of the both the modcl-bascd and subsample bootstrap fidelity as- 
sessment methods arc, to first order, given by the cost of the underlying calibration and imag- 
ing algorithm, scaled by the number of bootstrap resamples, A^^. We note above, however, 
that this algorithm is highly scalable in parallel computing environments as each resample 
can effectively be processed separately. 

In this numerical study, we have estimated imaging fidelity by exploring only the con- 
tribution of the thermal noise variation. A recent Bayesian approach to imaging fidelity 
assessment for optical astronomical images is presented by Esch et al. (2004). 

Bootstrap studies in general require limited statistical model assumptions and can be 
used directly with the measured data realization observed by a radio interferometer array. As 
such, they offer significant practical advantages, and wc believe them to be a viable approach 
to imaging fidelity assessment in the automated pipeline reduction of data from current and 
future radio interferometers. 



6. Conclusions 

On the basis of this numerical study, wc draw the initial conclusion that modern re- 
sampling techniques in computational statistics offer significant promise as imaging fidelity 
assessment techniques for common calibration and imaging processes used in radio interfer- 
ometry. This conclusion apphes both to the model-based bootstrap and subsample bootstrap 
methods. This computationally intensive approach is now tenable as a result of the recent 
geometric advances in available computing resources, and we believe these fidelity assessment 
techniques have an important role to play in automated pipeline reduction environments for 
modern radio telescopes. 

The polarization self-calibration algorithm developed as part of this study as the frame- 
work calibration and imaging test problem, offers a new approach to interferometric polar- 
ization self-calibration which makes no assumptions regarding the form of the polarization 
source model. 

Further work is required in assessing the statistical performance of the polarization 
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self-calibration algorithm and the resamphng techniques over a larger set of simulated data 
covering a broader range of observational parameters. 

This work was supported by a grant of computing time on the NCSA IA-32 Linux 
clusters under allocation AST 030025. We thank Drs R. Crutcher and T. Cornwell for 
reading the manuscript. We thank the referee for detailed comments on the manuscript 
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Table 1. Simulation D-term model values 



VLBA antenna 


\\D^\\ 


(deg) 




ct>^ (deg) 


Pie Town, NM 


0.03 


-57.0 


0.01 


23.0 


Kitt Peak, AZ 


0.02 


40.0 


0.01 


-5.0 


Los Alamos, NM 


0.03 


10.0 


0.05 


13.0 


Brewster, WA 


0.04 


87.0 


0.02 


3.0 


Fort Davis, TX 


0.09 


-15.0 


0.03 


56.0 


St. Croix, VI 


0.01 


-52.0 


0.01 


29.0 


North Liberty, lA 


0.05 


77.0 


0.01 


-23.0 


Owens Valley, CA 


0.03 


47.0 


0.01 


43.0 


Mauna Kea, HI 


0.03 


-7.0 


0.06 


83.0 


Hancock, NM 


0.03 


-27.0 


0.02 


35.0 



Note. — The D-terms tabulated here in a circulaj: basis translate 
to a D Jones matrix for antenna m of the form: 

\ j (22) 



Table 2. Simulation source model parameters 



Aa 


AS 






P.A. 




Q 


u 


V 


(mas) 


(mas) 


(mas) 


(mas) 


(deg) 


(Jy) 


(Jy) 


(Jy) 


(Jy) 


1 





0.2 


0.15 


30 


1 


0.2 


0.1 





2 -0.5 


-fO.l 


0.3 


0.1 


-10 


1 





0.25 


-0.1 



Table 3. Model-based bootstrap run parameters 



Run code 


max(JVj,) 


At(, (sec) 


Ml 


1 


60 


M2 


1 


150 


M3 


1 


300 


M4 


2 


900 




Note. — max(A'p) is the maximum 
degree of the model polynomial used in 
each bootstrap model interval Ati,. 
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Table 4. Subsample bootstrap run parameters 



Run code 


fs 


SI 


0.125 


S2 


0.25 


S3 


0.5 


S4 


0.75 



Note. — fe is 
the fraction of data 
deleted in each sub- 
sample realization 



Table 5. Bootstrap performance statistics 



Code 




'Vmse 
(xl0-i°) 


Ml 


1.00 


2.95 


M2 


1.01 


2.28 


M3 


1.00 


2.19 


M4 


1.02 


2.36 


Sl 


0.157 


6.43 


S2 


0.319 


5.15 


S3 


0.782 


3.42 


S4 


2.03 


3.04 


Note. — The bootstrap 
performance measures, Vf and 
vmse, are defined in the main 



text. 
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Table 6. MSE of CDF fit to parent distribution 



Code 


cmse 




(xlO-4) 


Ml 


11.1 


M2 


25.5 


M3 


5.35 


Mi 


S,()7 



Note. — The MSE 
of the CDF fit to 
the parent distribu- 
tion, Cmse, is de- 
fined in the main text. 
The value tabulated 
here is for visibility 
row 1000 in cross- 
correlation LR. For 
reference, the corre- 
sponding value from 
the Monte Carlo sam- 
ple is 3.48 X 10"*. 
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Fig. 1. — Simulation source model, convolved with a circular restoring 

beam of 156.007 fj,as, plotted in Stokes {/, Q, U, V}, using contour levels of 
{-64, -32, -16, -8, -4, -2, -1, 1, 2, 4, 8, 16, 32, 64} x 4.434 x 10"=^ Jy per 
beam. The peak brightness in Stokes {/, Q, U, V} is 4.434 x 10"^ 8.867 x 10"^ 9.591 x 10"^ 
and —3.836 x 10~^ Jy per beam respectively. 
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Fig. 2. — Parallactic angle variation at the Pie Town antenna for the data used in the Monte 
Carlo simulations. 
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Giga- wave length (U) 

Fig. 3. — itv-coverage for the data used in the Monte Carlo simulations. The simulated 
source has coordinates (q;2ooo = 12^^, ^2000 = +80'^'^^) 
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Self— call brotion Iteration number 



Fig. 4. — Mean-squared error in the off-diagonal D Jones matrix elements, solved for using 
polarization self-calibration, averaged over all antennas in the Monte Carlo simulation, and 
plotted on a logarithmic scale against self-calibration iteration number. The total calibration 
MSE is computed as: 2^ J2p£{R,L) Y^Zi^m " 'Ki){^m ~ ^m) ■> where is the instrumental 
polarization determined by the solver, and is the true value used to generate the simulated 
data. The instrumental polarization terms are defined in equation 7. 
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Fig. 5. — The measured imaging estimator bias, mse, and variance for Stokes Q, at polariza- 
tion self-calibration iteration numbers 1 and 10, for the Monte Carlo sample. Contour levels 
for the bias plots are at levels of {0.2n, n = —5, 5} of the local absolute maximum value 
(4.566 X 10~'^ at iteration 1, and 1.009 x 10~^ at iteration 10). Contour levels for the mse and 
variance nlots axe a,t iOAn. n= 1.10]- of the ]oca,] maximum. The nea.k mse is 2.155 x 10~^ 
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Fig. 6. — The measured imaging estimator bias, mse, and variance for Stokes U, at polariza- 
tion self-calibration iteration numbers 1 and 10, for the Monte Carlo sample. Contour levels 
for the bias plots are at levels of {0.2n, n = —5, 5} of the local absolute maximum value 
(4.476 X 10^'^ at iteration 1, and 8.261 x lO^'' at iteration 10). Contour levels for the mse and 

vflriannp nlnfs arp /fl In n = 1 1 fll- nf fVip Innal maYimnm T'Vip npnlr msp is 9 089 V 10"^ 
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Fig. 7. — The measured imaging estimator Stokes Q variance for the model-based bootstrap 
runs M1-M4, at the final polarization self- calibration iteration number 10. The contour 
levels are plotted at {O.ln, n = 1, 10} of the local maximum value, which are 1.666 x 10~^ 
(Ml), 1.852 X 10-6 (M2), 1.758 x 10"^ (M3), and 1.855 x 10"^ (M4) (Jy per beam)^. The 
parameters for the bootstrap run codes M1-M4 are defined in table 3. 
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Fig. 8. — The measured imaging estimator Stokes U variance for the model-based bootstrap 
runs M1-M4, at the final polarization self- calibration iteration number 10. The contour 
levels arc plotted at {O.ln, n = 1, 10} of the local maximum value, which are 1.737 x 10~^ 
(Ml), 1.765 X 10-6 (M2), 1.843 x 10"^ (M3), and 1.889 x 10"^ (M4) (Jy per beam)^. The 
parameters for the bootstrap run codes M1-M4 are defined in table 3. 
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Fig. 9. — The measured imaging estimator Stokes Q variance for the subsample bootstrap 
runs S1-S4, at the final polarization self-calibration iteration number 10. The contour levels 
are plotted at {O.ln, n = 1, 10} of the local maximum value, which are 1.827 x 10~^ (SI), 
3.849 X 10"^ (S2), 7.994 x 10"^ (S3), and 1.890 x 10"^ (S4) (Jy per beam)^. The parameters 
for the bootstrap run codes S1-S4 are defined in table 4. 
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Fig. 10. — The measured imaging estimator Stokes U variance for the subsample bootstrap 
runs S1-S4, at the final polarization self-calibration iteration number 10. The contour levels 
are plotted at {0.1?7., n = 1,10} of the local maximum value, which are 1.817 x 10~^ (SI), 
3.737 X 10-^ (S2), 8.262 x 10"^ (S3), and 2.100 x 10"^ (S4) (Jy per beam)2.The parameters 
for the bootstrap run codes S1-S4 are defined in table 4. 
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Fig. 11. — The variance scaling factor, Vf, plotted against the subsampling delete fraction, 




Fig. 12. — The empirical CDF for visibility row 1000, polarization cross-correlation LR, for 

the model-based bootstrap resamples Ml to M4. The CDF for the real part of the visibility 
data is plotted in the left column, and the CDF for the imaginary part of the visibility data 
in right column. The solid line is the expected normal CDF from the parent distribution. 
All data have been centered by subtracting the expected visibility value Vmn at this visibility 
row number . 
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Fig. 13. — The empirical CDF for visibility row 1000, polarization cross-correlation Li?, for 
the Monte Carlo sample. The CDF for the real part of the visibility data is plotted in the 
left column, and the CDF for the imaginary part of the visibility data in the right column. 
The solid line is the expected normal CDF from the parent distribution. All data have been 
centered by subtracting the expected visibility value Vmn at this visibility row number. 



